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ABSTRACT 

We develop a formalism for studying the dynamics of massive black hole binaries embedded in 
gravitationally-bound stellar cusps, and study the binary orbital decay by three-body interactions, the 
impact of stellar slingshots on the density profile of the inner cusp, and the properties of the ejected hyper- 
velocity stars (HVSs) . We find that the scattering of bound stars shrinks the binary orbit and increases 
its eccentricity more effectively than that of unbound ambient stars. Binaries with initial eccentricities 
e ^ 0.3 and/or unequal-mass companions (M 2 /Mi <, 0.1) can decay by three-body interactions to the 
gravitational wave emission regime in less than a Hubble time. The stellar cusp is significantly eroded, 
and cores as shallow as p oc r~ may develop from a pre-existing singular isothermal density profile. A 
population of HVSs is ejected in the host galaxy halo, with a total mass ~ M2. We scale our results to 
the scattering of stars bound to Sgr A*, the massive black hole in the Galactic Center, by an inspiraling 
companion of intermediate mass. Depending on binary mass ratio, eccentricity, and initial slope of the 
stellar cusp, a core of radius ~ 0.1 pc typically forms in 1-10 Myr. On this timescale about 500-2500 
HVSs are expelled with speeds sufficiently large to escape the gravitational potential of the Milky Way. 

Subject headings: black hole physics - methods: numerical - stellar dynamics 



1. INTRODUCTION 

In the standard paradigm of cosmic structure formation, 
it is expected that many wide massive black hole binaries 
(MBHBs) will form following the mergers of two massive 
galaxies (e.g. Begelman, Blandford, & Rees 1980; Volon- 
teri, Haardt, & Madau 2003; Mayer et al. 2007). The 
binary will subsequently shrink due to stellar or gas dy- 
namical processes and may ultimately coalesce by emitting 
a burst of gravitational waves. It was first proposed by 
Ebisuzaki, Makino, & Okumura (1991) that the heating 
of the surrounding stars by a decaying SMBH pair would 
create a low-density core out of a preexisting cuspy (e.g. 
p cx r -2 ) stellar profile. In a purely stellar background a 
'hard' binary shrinks by capturing the stars that pass close 
to the holes and ejecting them at much higher velocities, a 
super-elastic scattering process ('gravitational slingshot') 
that depletes the nuclear region. Observationally, there is 
evidence in early-type galaxies for a systematically differ- 
ent distribution of surface brightness profiles, with faint 
ellipticals showing steep power-law profiles (cusps), while 
bright ellipticals have much shallower stellar cores (e.g. 
Faber et al. 1997; Ravindranath et al. 2001). Dwarf el- 
lipticals seem to elude this paradigm showing somewhat 
flat profiles, similarly to bright ones (Graham & Guzman 
2003). Late-type spirals tend to show steep central cusps, 
as in the case of the Milky Way or Andromeda. Detailed 
N-body simulations have confirmed the cusp-disruption ef- 
fect of a hardening MBHB (Makino & Ebisuzaki 1996; 
Quinlan & Hernquist 1997; Milosavljevic & Merritt 2001), 
while semi-analytic modelling in the framework of hier- 
archical structure formation theories has shown that the 
cumulative damage done to a stellar cusps by decaying 
black hole pairs may explain the observed correlation be- 



tween the 'mass deficit' (the mass needed to bring a flat 
inner density profile to a r~ 2 cusp) and the mass of the 
nuclear black hole (Volonteri, Madau, & Haardt 2003). 

This is the third paper of a series aimed at the detailed 
study of the interaction between MBHBs and their stellar 
environment. In Sesana, Haardt, & Madau (2006,2007a, 
hereafter Paper I and Paper II), we analyzed the three- 
body scatterings between a MBHB and background stars 
unbound to the binary. The assumption of a fixed back- 
ground breaks down once the binary has ejected most of 
the stars on intersecting orbits, and the extraction of en- 
ergy and angular momentum from the binary can continue 
only if new stars can diffuse into low-angular momentum 
orbits (refilling the binary's phase-space "loss cone"), or 
via gas processes. In galaxies with inner cores or shal- 
low cusps, only a small fraction of the loss cone is con- 
fined within the sphere of influence of the binary, and the 
approximation of a background of unbound stars is rea- 
sonable. A similar argument holds also in the case of 
a galaxy with a steep cusp hosting a nearly equal-mass 
binary (Mi ~ M2). The radius of influence of such a 
pair, fjnf = G(Mi + M2)/(2a 2 ) where a is the stellar ve- 
locity dispersion, is of the order of the binary hardening 
radius, = GM2/ '4<r 2 , and only few low-angular momen- 
tum stars have orbits with semi-major axis <J r; n f . This is 
not true for unequal-mass binaries, where r; n f 3> dfe, and 
almost all interacting stars are bound to Mi. 

In this paper we develop a formalism for studying the 
dynamics of MBHBs embedded in gravitationally-bound 
stellar cusps. The plan is as follows. In § 2 we describe 
our suite of three-body scattering experiments between the 
black hole pair and ambient bound stars. Our numerical 
results are used in § 3 and § 4 to construct a "hybrid 
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model" of binary dynamics and investigate the orbital de- 
cay and shrinking of MBHBs in time-evolving stellar cusps. 
The properties of the ejected HVSs are discussed in § 5. 
The massive black hole (Sgr A*) in the Galactic Center 
and the stars around it offer a unique opportunity to study 
stellar dynamics in the extreme environment around a rel- 
ativistic potential. The scattering of stars bound to Sgr 
A* by an inspiraling intermediate-mass black hole (IMBH) 
is treated in § 6. Finally, we present a brief summary in 
§7. 

2. SCATTERING EXPERIMENTS WITH BOUND STARS 

Consider a MBH of mass Mi surrounded by a stellar 
cusp (with density profile p oc r~ 7 , 7 > 0), interacting 
with a secondary "intruder" hole of mass M2 < M\ . As 
the binary separation decays, the effectiveness of dynam- 
ical friction slowly declines because distant stars perturb 
the binary's center of mass but not its semi-major axis. 
The bound pair then loses orbital energy by capturing am- 
bient stars and ejecting them at much higher velocities, a 
three-body scattering process known as the 'gravitational 
slingshot'. For unequal-mass binaries (mass ratio q <; 0.1), 
o-h *C ^jnf, implying that the stellar mass inside is 
M cusp (< ah) -C Mi. In this case, the contribution to the 
potential energy given by the stellar distribution during a 
binary-star interaction can be ignored, and the problem 
can be tackled by means of three-body scattering experi- 
ments. Following the guidelines fully described in Paper I, 
we performed numerical experiments to study the interac- 
tion of a MBHB with a star of mass to* (to* <C M% < Mi) 
bound to Mi. Note that, for extreme mass ratios q <C 1, 
the encounter is essentially a two-body scattering as the 
star and the secondary hole move in the static potential of 
the primary. 

2.1. Initial conditions and orbit integration 

The integration of the three-body encounter equations 
is greatly simplified by setting the center of mass of the 
binary at rest in the origin of the coordinate system. The 
binary orbits counterclockwise in the [x, y) plane, and the 
apoastron of M2 is located along the positive x-axis. Stars 
are drawn from a spherical isotropic distribution bound to 
Mi. The problem is completely defined by ten variables: 

1. the binary mass ratio q = M2/M1; 

2. the binary eccentricity e; 

3. the stellar mass to*; 

4. the initial distance between the primary hole and the 
star, r* = |r*|; 

5. the specific energy of the stellar orbit around Mi, 
E* = — GMi/(2o*), or equivalently the semi-major 
axis of the stellar orbit, a*; 

6. the specific angular momentum of the stellar orbit 
L„ = |v x r*|, where v = |v| is the stellar velocity 
relative to Mi; 

7. four angles: 9 and (f> describing the latitude and lon- 
gitude of the star, tjj defining the orientation of v± 
(the stellar velocity component normal to r*) i.e. se- 
lecting the orbital plane of the star, and \P the initial 
binary phase. 
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Fig. 1. — Geometry of the three-body scattering experiments. 
The stellar orbits is determined by E* and L* , while its orientation 
with respect to the binary plane is determined by angle tp. See text 
for definition of all symbols. 

Compared to the unbound scattering problem (which 
is defined by nine variables, see Paper I), there is here 
the extra variable r*, and and L* replace the asymp- 
totic initial speed and impact parameter of the incoming 
field star. A sketch of the experiment set-up is given in 
Figure 1. In each numerical integration, the binary mass 
M = Mi + M2 and separation a are set equal to unity, 
and the stellar mass to m* — 10™ 6 Af. Scattering events 
are simulated in sets of 5 x 10 4 trials for fixed values of 
q and e. We sample five values of the binary eccentric- 
ity, e = 0,0.1,0.3,0.6,0.9, and five values of the binary 
mass ratio, q = 1/9, 1/27, 1/81, 1/243, 1/729, for a grand 
total of 25 models. The distribution of stars bound to 
Mi is modeled as follows. The stellar semi-major axis a* 
is randomly sampled from fifty logarithmic bins spanning 
the range 0.03a < a* < 10a. The angular momentum is 
sampled in the interval [0, L* max ] according to an equal 
probability distribution in L%, where L* max = GMia* is 
the specific angular momentum of a circular orbit of ra- 
dius a*. A population of stars with such distribution in 
L* has mean eccentricity (e*) = 0.66, corresponding to 
(v±) = 2(v||). This condition defines an isotropic stellar 
distribution (e.g. Quinlan, Hcrnquist, & Sigurdsson 1995). 

The quantities E* and L* define the shape of the stellar 
orbit. We sample the initial value of r* from the distribu- 
tion P(r*)dr*, 
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which is proportional to the the fraction of time spent by 
the star at distance between r* and r* + dr* from Mi . The 
probability "P(r*) is defined in the range r*__ < r* < r* )+ , 
where 



r*,_ = L 2 J(GM 1 + y ggg + 2EM ) 
r*,+ = L 2 J(GMi - y/G 2 Mf + 2E*Ll). 



(2) 



The angles 9 and <fi are randomly generated to reproduce a 
uniform density distribution over a spherical surface cen- 
tered on Mi , while the orientation angle ip is chosen from 
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a uniform distribution in the range [0, 2tt]. We start nu- 
merical integration with M 2 at its apoastron (5* = 0), and 
have checked that the chosen initial phase of the binary 
does not affect our results. The orbit of the pair is in the 
x — y plane with the center of mass at coordinates (0, 0, 0). 

The nine coupled, second order, differential equations 
of the three-body problem are integrated using the most 
recent version of the subroutine DOPRI5, based on an ex- 
plicit Runge-Kutta method of order 4(5) due to Dormand 
& Prince (1978). A complete description of the integrator 
can be found in Hairer, Norsett, & Wanner (1993). The 
integration is stopped if any of these events occurs: 



the star leaves the sphere of radius r = 
[10 10 /i/M] 1 / 4 ao with positive total energy, where 
[i = M1M2/M is the reduced mass of the binary 
and clq its initial semi-major axis. At r the force 
induced by the quadrupole moment of the binary is 
10 orders of magnitude smaller than the total force 
acting on the star at a distance ao; 



the physical integration time exceeds 1 Gyr; 
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Fig. 3. — Example of a weak three-body interaction. Left panel: 
projection onto the binary (x — y) plane. Right panel: projection 
onto the y — z plane 



• the integration reaches 10 8 time steps. This typically 
corresponds to the complete integration of 10 4 — 10 5 
binary orbits. 




Fig. 2. — Example of a strong three-body interaction unbind- 
ing the star. The parameters of the system are listed in the fig- 
ure. The star is ejected towards the lower right with a kick velocity 
V ~ 0.05V C . The thick circle marks the orbit of M2, while the small 
central dot marks the location of M\. The insets zoom-in at the 
moment of the scattering between the star (approaching from the 
left) and M2 (approaching from below). 



Examples of integrated stellar orbits are given in Fig- 
ures 2 and 3. Figure 2 shows a strong interaction: the 
star has a close encounter with M2 and leaves the binary 
sphere of influence with positive energy. Depending on 
the value of a» we have a variable number of numerically 
"abandoned events" , i.e., bound triple systems that do not 
readily produce an ejection. Figure 3 shows one of such 
events: the periastron of the star is much larger than the 
binary orbital separation a, and the star is not perturbed 
by the binary quadrupole field, except for the precession 
of its periastron. 



2.2. Tests 

We have performed a number of tests to check the sen- 
sitivity of our results on numerics. Because of the in- 
trinsically chaotic nature of the three-body problem, the 
properties of the ejected stars are meaningful only for a 
statistically-significant sample. The integration of the full 
three-body problem allows us to directly control the con- 
servation of total energy and angular momentum. The 
code adjusts the integration stepsize to keep the frac- 
tional error per step in position and velocity, e, below 
10 -13 . This allows a total energy conservation accuracy 
of AE/E ~ 10~ 9 in a single orbit integration, i.e. for 
m*/M ~ 10~ 6 the energy of the star is conserved at level 
of one part in a thousand during a single orbit. We varied 
e between 10~ n and 10~ 15 and m*/M between 10~ 5 and 
10 -7 , and found no significant differences in the statistics 
of the ejected population. We have also checked that the 
mean energy and angular momentum exchanges scale lin- 
early with m*. A set of longer orbit integrations was per- 
formed to test that the number of abandoned events does 
not depend on the finite numbers of timesteps allowed, 
again finding no systematic decrease in such a number. 
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2.3. Outputs 

Each logarithmic bin in a* was sampled by 10 3 stars. 
We calculated the mean energy and angular momentum 
exchange between the MBHB and the stars, the final ve- 
locity and angular distribution of the scattered stars, the 
fraction / e j of stars that are ejected in the interaction, and 
the ejection timescales. The fraction / e j is plotted versus 
a* I a in Figure 4 for different values of the binary param- 
eters q and e. Eccentric binaries can eject stars that are 
initially very tightly bound to Mi, i.e. with tts|s CIS small 
as 0.1a. For a* ^ 3a, the expelled fraction declines signif- 
icantly, regardless of binary eccentricity. It is also evident 
that the ejection process is more efficient for large values 
of q. 
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Fig. 4. — Fraction of stars ejected in the interaction as a function 
of a t /a,. Thin curves are for q = 1/9 and different eccentricities: 
e = (solid line), e = 0.1 (long-dashed line), e = 0.3 (short-dashed 
line), e = 0.6 (dot-dashed line), e = 0.9 (dotted line). The thick 
dashed line is for q = 1/729 and e = 0.3. 

Figure 5 shows the mean fractional eccentricity change 
of the pair after each scattering, (Ae/e), as a function of 
a*/a. This quantity is found to scale linearly with m*/M. 
Stars with a» < a typically tend to reduce the binary ec- 
centricity, while stars with a* > a work in the opposite 
direction. An eccentric binary spends most of its period 
near its apocenter, so in the case a* > a the probabil- 
ity of a close star-binary encounter (and subsequent star 
ejection) is maximal at binary apocenter. In the instanta- 
neous interaction the binary velocity decreases. As close 
to the apocenter the binary velocity is nearly tangential, 
the binary is forced on a more radial orbit. On the con- 
trary, a star with a» < a "feels" the secondary hole M 2 
when this approaches the primary Mi at a distance ~ a* . 
At that point, the interaction with the star is unlikely 
to occur close to the pericenter of M 2 (being the time 
spent by M 2 at the pericenter very small), and typically 
also extracts a large radial component from the velocity 



of the secondary black hole, hence causing circularization. 
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Fig. 5. — Mean fractional eccentricity change, (Ae/e), of the bi- 
nary after a scattering, as a function of a* /a. The assumed initial 
eccentricity is 0.6, and the different curves are for q = 1/729 (solid 
line), q = 1/243 (long-dashed line), q = 1/81 (short-dashed line), 
q = 1/27 (dot-dahed line), and q = 1/9 (dotted line). Note the 
different scales of the positive and negative j/-axis. 
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Fig. 6. — Differcntail distribution of ejection timescales r. The 
assumed initial eccentricity is 0.3, and the different curves, from top 
to bottom, are for q = 1/9, 1/27, 1/81, 1/243, and 1/729. 
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Note that in this case the absolute value of Ae is larger, 
according to the aj a* dependence in equation (9) (see the 
different positive and negative j/-axis scales in Fig. 5). We 
will return on this point in Section 3, giving a simple, phys- 
ically motivated, mathematical derivation of this qualita- 
tive argument. 

Figure 6 depicts the fractional number of stars ejected 
in the time interval r and r + dr as a function of the 
ejection timcscale r. The latter is measured in units of 
the binary orbital period at separation ao = a(t = 0), 
Po = 2TTy/a^/(GM), and is defined as the time elapsed 
from the start of numerical integration to the moment the 
interacting star reaches, with positive energy, a distance 
3> a from the binary center-of-mass. For large values of q, 
it typically takes the MBHB just a few orbits to expel the 
star, while lowering q makes the distribution of slingshot 
timescalcs broader and flatter. The ejection rate remains 
approximately constant, or decreases slowly, for r < 5/g 2 , 
and drops dramatically afterwards. The binary eccentric- 
ity plays no role in all the tested cases. Figure 7 shows the 
mean ejection timescale (i.e., the integral of T(df e j/dr) in 
dr) as a function of the star eccentricity e*, for different 
values of q and e. In first approximation, one can assume 
that the cusp remain isotropic during the binary orbital 
evolution. 
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Fig. 7. — Mean ejection timescale < r > of stars as a function of 
the star eccentricity e», assuming a MBHB with e = 0.1 (solid line) 
and e = 0.9 (dashed lines). The two upper (lower) curves are for 
g = l/81 (g = l/9). 

We checked that our results are practically the same con- 
sidering 5a rather than 10a as the outer boundary in the 
sampled a* distribution. This is because the binary-stars 
interaction is dominated, in terms of energy and angular 
momentum exchange, by those stars whose a* is <~ a. 

From our scattering experiments, we can finally de- 
rive the bivariate distribution functions hi(V, 9\a*) and 
hi(V, </>|a»), along the same line described in Paper II. For 



a given a*, we record the number of of stars with ejec- 
tion speed in the interval V, V + dV , leaving the binary 
with latitude (longitude) at infinity in the interval 9,9 + d9 
(4>, (j) + d(f>) . The differential distributions are normalized 
as follows: 

i>OC i>7T />OC pIlT 

/ dV d9h 1 (V,9\a*) = dV #/n(V, <£K) = 
Jo Jo Jo Jo 

(3) 

The subscript "1" is meant to indicate that the scatter- 
ing experiments are performed for a binary at separation 
a = l. 

3. HYBRID MODEL OF BINARY DYNAMICS 

To study the impact of the gravitational slingshot on 
the dynamical evolution of a MBHB, we have developed 
a self-consistent hybrid model, in which numerical results 
of scattering experiments in a fixed stellar background are 
coupled to an analytical formulation of loss-cone depletion. 
This technique allows us to simultaneously follow the or- 
bital decay of the pair as well as the time evolution of 
the stellar cusp. The hybrid model is similar, in spirit, to 
the case extensively discussed in Paper II of an unbound 
stellar population interacting with a black hole pair. Scat- 
tering experiments are performed for a binary at fixed or- 
bital separation a, and the results scaled to any value of a. 
Then, by specifying the pace at which a evolves, we can 
solve for the binary orbital decay with time. Mathematical 
details, however, are different from the scheme developed 
in Paper II. In the unbound case, the pair interacts with 
a stellar population approaching its sphere of influence at 
a given rate, and the typical three-body interaction is a 
"fast" process. For stars bound to Mi instead, the inter- 
acting stellar population is in place from the beginning. 
The temporal evolution of the system is then determined 
by the ejection timescale rather than by supply rate. 

It is convenient to describe the stellar cusp using 
dN^/da*, the differential number of stars orbiting Mi with 
semi-major axis in the range a*, a* + da». For an isotropic 
stellar distribution 



dN* 
da* 



C 4:iralp(r = a»), 



where the stellar density profile is 



P(r) = Po 



(4) 



(5) 



and C is a fudge factor that depends on the cusp slope, 
C = (1.013,0.872,0.831) for 7 = (2,1.75,1.5) (Ivanov, 
Polnarev, & Saha 2005). Two differential equations de- 
termine the rate of change of orbital separation and ec- 
centricity: 



da 
~dt 



2a 2 



GM1M2 



A£- -4 da*, 

aa*ar 



and 



de 
~dt 



Ae- — -^da*. 
da*dt 



(6) 



(7) 



We start numerical integration at time t — 0, orbital sep- 
aration ao, and binary eccentricity eo- The terms Ae, AS, 
and d 2 N e j/da*dt are measured from our scattering exper- 
iments. 
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The scaling of the eccentricity change with binary pa- 
rameters m*/M, e, q, and a/a* can be derived from (Quin- 
lan 1996) 



Ae = - 



(1 



2e 



AS 
~S 



2AC Z 



(8) 



where C z = \i^J GMa(l - e 2 ) and £ = -GM 1 M 2 /(2a) are 
the total angular momentum and energy of the binary, re- 
spectively. From AS = — AS* ~ — GM\m*j(2a*) (in a 
typical encounter, on average, the star gets a kick ~ V c ), 
and AC Z = —AC* ~ — C* = -m.y'GMfl.fl — e 2 ), as- 
suming Mi ~ M, we have 



a M 
Ae — 

m* 




(9) 



This conhrms the qualitative argument developed in Sec- 
tion 2.3 and accounts for the features shown in Figure 5: 
the mean Ae scales with the inverse of the mass ratio, and, 
from < e* >~ 0.67, results negative for a*/a < (1 — e 2 ) 1 / 3 
(promoting circularization), and positive for larger a* /a 
(hence increasing the binary eccentricity). 

The binary energy change per scattering, AS, is related 
to the change of stellar energy AS* by AS* ~ — AS, where 
we have neglected the energy change of the binary center 
of mass (which is a factor m*/M smaller than AS). The 
quantity AS* can be written as 

AS*(a*,a, e) = GM ^ m * + Ifc m *V 2 (a), (10) 



where V c (a) = \J GM/a is the circular velocity of the bi- 
nary. The numerical factor k depends upon the ratio a* /a, 
and e, and is derived from our numerical experiments. 

The term d 2 N e -Jda*dt quantifies the number of stars 
orbiting M\ within a* and a* + da* that are ejected from 
the system in the time interval between t and t + dt. This 
term depends on a*, a, e and t, and is determined as fol- 
lows. From scattering experiments we derive the distri- 
bution P(t, s\e)dr describing the probability that a star 
at s = a* /a becomes unbound from a binary with eccen- 
tricity e in the time interval (t, t + dr). Note that the 
distribution of ejection timescales plotted in Figure 6 is 
simply the probability function P averaged over s, 



or 



/ P(r, s\e)ds 

7^ ' 



(ii) 



When s<Cl (the exact value depending on q and e) most 
stars remain bound to M\, and P = (see Fig. 4). As 
s — > 1 all interacting stars are instead expelled. We can 
see from Figure 6 that, in the ejection regime, the dis- 
tribution P is nearly independent on s and has the same 
functional form (as a function of r) of df e j/dr. From P we 
can then compute d 2 N e - 3 /da*dt as a function of t if we set 
t = (a /a) 3 / 2 t, where the term (a /a) 3 / 2 accounts for the 
change in the time units of P as a changes. Though not 
formally correct, the scheme catches the basic physics of 
the interaction in a time-dependent fashion. 

Suppose now that in a small time interval Ato after the 
beginning of the interaction the binary remains at constant 



separation ao corresponding to so = a*/a>o. The number 
of stars (with semi- major axis in the interval a*, a* + da*) 
ejected in such time interval is 



d2Nci :(0)At = ^P( , So )At , 



da*dt 



da* 



(12) 



where for simplicity we have omitted the dependence of 
P on the eccentricity. After Ato and for an interval Ati, 
the binary settles to a new separation ai, corresponding 
to si = a*/a\. The number of stars ejected in the time 
interval Ato, Ato + At\ is 



dN* 
da* 



-To 



P(At Q , Sl ) (ao 

SZ o nt',s,)dt' 'U 



3/2 



The integral at the denominator of the right-hand side 
renormalizes the distribution P in the time interval 
[Ato, oo ] so that the correct number of stars is involved. 
Iterating we have 



dN* 
da* 



where t = Yli=o ^*»- l n differential form: 



p(t, Sj ) 



At, 



i=0 



f~P(t>, Sj )dt' J \a 



ao 



d 2 N cj 
da*dt 



(t) 



dN* 
da* 



f d 2 N ci 
Jo da*dt' 



dt' 



P(t,s) 



3/2 



(14) 



ao 



f t °°P(t>,s)dt> [a(t) 
(15) 

where s = s(t) = a*/a(t). Simultaneous numerical inte- 
gration of the tree coupled equations (6), (7) and (15) self- 
consistently solves for the evolution of the binary and the 
depletion of the stellar cusp. The integration is performed 
using the subroutine DOPRI5. The time-step is adapted 
in order to keep the fractional error per step < 10~ 10 , a 
value much lower than the error associated to the employed 
linear interpolation of Ae and AS. 

The bivariate distributions h\{V,0\a*) and hi(V,4>\a*) 
derived from our suite of scattering experiments (see 
§ 2) can be convolved with the ejection rate d 2 N e j/da*dt 
to compute the final velocity distributions h(V, 6) and 
h(V, (j)). The procedure is similar to that described in the 
Appendix of Paper II. As the binary shrinks to separation 
a < 1, the normalized distribution of ejection velocities for 
stars with semi-major axis in the interval a*, a* + da* is 



(16) 



h a (V,e\a*) = 4=^1 

where the ejection speed V was shifted by the factor 1 / ^fa 
to account for the increase in the circular velocity of the 
binary V c as the pair shrinks. The prefactor l/\/a normal- 
izes the distribution according to equation (3). The kick 
velocity function of the expelled population as a whole can 
then be written as 



L af So 00 ^a(V,e\a*)da*da 



h(V,6) = 



af poc p-rr poo d 2 N c : 



rdf POO (* 7T (* C 

Ji Jo Jo Jo 



da*da 



h a {y, 9\a*)da*dad0dV 



where a/ is the final binary separation, 

d 2 N ci _ d 2 N e j dt 
da*da da*dt da' 



(17) 



(18) 
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and dt I da is given by equation (6). The distribution 
h\(V, 6\a*) is evaluated for different values of the eccen- 
tricity e. We can account then for the evolution of binary 
eccentricity by interpolating the hi distribution on a grid 
of e-values as the orbit decays. The above procedure re- 
turns the velocity distribution in units of V c (ao). The cal- 
culation of h(V, (/)) can be performed following the same 
lines. 

Finally, integrating the rate d 2 N e -Jda*dt over the en- 
tire evolutionary history of the binary yields the quantity 
dNcj/da*, which can be compared to dN/da* to study the 
depletion of the stellar cusp by three-body interactions. 
The Monte Carlo technique developed to address this point 
is described in § 5.1. The functions dN e -Jda* and h(V, 6) 
can be used to check energy conservation, by simply equat- 
ing the total energy gained by the stars to the increment 
in the binding energy of the black hole pair, i.e. 



m, 



f 

Jo 



oo dJV ei 



I 1 f c 
+ 2 TO * Jo 



$(a») da* + 
V 2 N ci h(V) dV 



GM t M 2 
2 



\ a f a o j ' 



(19) 

where $(a„) = GMi/a*, h(V) = J h(V,0)dO, and N ci = 
J(dN C j/da*) da*. We have checked that energy is con- 
served to better than <~ 1%. 

It is important to remark at this point that equations 
(6) and (7) are independent of the absolute value of M and 
m*. Indeed, in equation (6), the total number of interact- 
ing stars is oc M 2 /m >t , and the energy exchange term is 
oc Mm*. The scaling factor I/M1M2 cancels out, so that 
the orbital evolution depends only on q and e. The same 
consideration holds for equation (7): as Ae oc to*/M while 
the number of interacting stars oc M 2 / to* , the dependence 
on M and m* cancels out. 

4. BINARY EVOLUTION 

For comparison with our previous results on the scatter- 
ing of unbound stars (Paper II), it is convenient to intro- 
duce the stellar velocity dispersion a, and model the outer 
stellar component as a singular isothermal sphere (SIS) 
with density profile 



p(r) 



(20) 



2irGr 2 ' 

We assume that this profile extends inward up to the char- 
acteristic radius r within which the total stellar mass is 
2 Mi (Merritt 2004). Matching the inner cusp described 
by equation (5) to the outer SIS at r yields 

(3-7)^. (2D 



ro 



a' 



We also assume that stellar-binary interactions start at 
separation ag where the enclosed mass stellar is M*(< 
a ) = 2M 2 , yielding 

1/(3-7) 



a 



1 + 9 



(3-7) 



GM 



(22) 



This assumption is motivated by recent N-body simula- 
tions of the hardening of unequal MBHB in stellar cusps. 
Matsubayashi, Makino, & Ebisuzaki (2007) found that dy- 
namical friction is efficient in driving orbital decay only as 
long as the stellar mass inside the binary semi-major axis 
is > M 2 . Beyond this point, the evolution of the pair is 
driven by three-body interactions with individual stars. 

4.1. Orbital evolution and eccentricity growth 

Our results on binary orbital decay are summarized in 
Table 1. A MBHB can shrink by factors ranging from 6 
to 18 depending on q, eo, and 7. The decay factor ao/o/ 
grows with eo, as more eccentric binaries can expel stars 
that are more tightly bound (see Fig. 4). The importance 
of the initial eccentricity eo is modest for low values of 7, 
since in shallow cusps the number of tightly bound stars is 
small anyway. The factor ao/a/ is also a weakly increasing 
function of q, as for higher q the radius of influence of the 
binary is larger (in terms of a*/a, see Fig. 4). The de- 
pendence of the decay factor on 7 simply reflects the fact 
that, for steeper stellar cusps, the mean binding energy of 
the stars is larger. 4 
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1.5 


7 = 


1.75 


7 = 


2 


q 


eo 


ao/a,f 


e.f 


ao/a/ 


e.f 


ao/a/ 


e f 




0.1 


9.63 


0.608 


9.89 


0.350 


11.38 


0.179 


1 

9 


0.5 


9.99 


0.972 


11.55 


0.907 


15.77 


0.753 


0.9 


10.06 


0.998 


11.80 


0.992 


17.73 


0.969 




0.1 


8.06 


0.691 


8.35 


0.532 


10.19 


0.408 


1 

27 


0.5 


8.26 


0.959 


9.35 


0.862 


12.35 


0.710 


0.9 


8.27 


0.996 


9.64 


0.988 


14.03 


0.958 




0.1 


6.99 


0.755 


7.75 


0.650 


9.39 


0.542 


1 

81 


0.5 


6.90 


0.922 


7.81 


0.828 


10.14 


0.717 


0.9 


6.89 


0.996 


7.81 


0.974 


11.00 


0.937 




0.1 


6.49 


0.906 


7.28 


0.805 


9.30 


0.688 


1 

243 


0.5 


6.39 


0.971 


7.25 


0.914 


9.92 


0.818 


0.9 


6.38 


0.962 


7.19 


0.986 


10.09 


0.955 




0.1 


6.12 


0.881 


6.91 


0.814 


8.94 


0.724 


1 

729 


0.5 


5.95 


0.919 


6.91 


0.869 


9.09 


0.797 


0.9 


5.94 


0.977 


6.92 


0.953 


9.41 


0.900 



Table 1 

Binary shrinking factors and final eccentricities from the hybrid 
model. 

Binary eccentricity also grows as a function of q and 7. 
A shallow cusp increases the relative importance of stars 
with large a*, and the eccentricity growth is then larger. 
Moreover, binaries with large q are more effective in eject- 
ing stars with (1,/aCl. As these stars act to circularize 
the binary orbit, nearly equal-mass binaries decay follow- 
ing less eccentric orbits. Figure 8 shows examples of binary 
evolution as a function of time. Note the different scale of 



4 Note that we have not included the effect of "returning stars" on the decay of the binary. These are ejected stars that do not escape the 
host bulge, return on small impact parameter orbits, and can have a secondary super-clastic interaction with the MBHB. The role of secondary 
slingshots was analyzed in Paper II, where it was found that they can boost orbital decay by as much as a factor of 2 for nearly equal-mass 
binaries, but do not contribute significantly (less than 20% for q = 1/9 and less than 5% for q = 1/81, see Table 1 of Paper II) to binary 
hardening for very unequal-mass pairs. In the bound case, the effect of returning stars is likely to be even smaller. Consider for example the 
case of an isothermal stellar profile: the impact of secondary slingshots is proportional to the number of stars that can interact with the MBHB 
more than once, i.e. to the size of the loss-cone after the first interaction. For bound stars this is at least a factor of ~ 2 smaller than in the 
unbound case: so even for q = 1/9, returning stars would cause at most a ~ 10% increase in binary hardening. 
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the time axis for the three mass ratios: as shown in Fig- 
ure 6, it takes nearly equal-mass pairs a smaller number 
of binary orbits to unbind the stellar cusp. 

Results for the 7 = 2 cusp can also be compared to 
those obtained in Paper II for the case of a binary in- 
teracting with a population of unbound stars. (In Pa- 
per II, shrinking factors were normalized to the "harden- 
ing radius" a/j = GM 2 /4er 2 , where, by definition, a = 
4a?j for 7 = 2.) Binaries with e = 0.1 and q = 
(1/9, 1/27, 1/81, 1/243), for example, shrink by the factors 
Oft/a/ = (2.85, 2.55, 2.35, 2.33) according to Table 1, com- 
pared to the corresponding ah/df — (2.09, 1.49, 1.19, 1.06) 
for the unbound case, i.e., very unequal-mass binaries can 
only decay by extracting the cusp binding energy. 
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Fig. 8. — Time evolution of the binary semi-major axis (upper 
panels) and eccentricity (lower panels). Left panels: q = 1/9. Cen- 
tral panels: q = 1/81. Right panels: q = 1/729. In all panels the 
thick lines are for a 7 = 2 stellar cusp, while the thin lines are for 
7 = 1.5. The initial binary eccentricity is assumed to be 0.1 (solid 
lines) and 0.5 (dashed lines). 

The evolution of the eccentricity in the bound and un- 
bound cases is compared in Figure 9, for two different 
initial values of e . Compared to the scattering of un- 
bound stars, the binary eccentricity grows to larger values. 
This can be understood by the following argument. In the 
bound case treated here, there is a great deal of stars able 
to extract angular momentum from the binary right from 
the start. Then, for large a, e shows a steep increase. 
When the orbit shrinks, however, stars with a* < a (in- 
creasing circularization rather than eccentricity, see Fig. 5) 
become more and more relevant for the binary evolution as 
stars with a* > a get progressively depleted, and, as result, 
e increases at a reduced pace. In the unbound case, the 
binary interacts with the same distribution of stars inde- 
pendently upon its separation. As, in general, significant 
changes of e happen whenever a strong star-binary interac- 
tion occurs, the eccentricity grows more rapidly for small 



a because stars are slower (in terms of binary circular ve- 
locity), and hence more easily captured in "quasi" bound 
orbits before ejection. As a matter of fact, the eccentric- 
ity growth rate K (see Paper I) increases for a <C a/,. In 
this separation regime, most of the stars interacting with 
the binary are captured in temporarily bound orbits, and 
would experience a dynamical interaction similar to that 
of bound stars. 



<u 




1 0.4 

a / a h 

Fig. 9. — Evolution of binary eccentricity as a function of orbital 
separation, for the case of interactions with unbound stars (Paper I 
II, dashed lines), and with a bound cusp (solid lines). The binary 
is embedded in an SIS, its mass ratio is q = 1/9, and its initial ec- 
centricity is eo = 0.3 and 0.6. Note that the scattering experiments 
with unbound stars start at a = a^, while those with bound stars 
start at a = Aa^ (see text for details). In the unbound case the 
evolution of eccentricity is negligible (K Ri 0, Scsana ct al. 2006) 
for a > a^- 

4.2. Final coalescence 

Consider now, as in Paper II, the separation at which 
a MBHB can coalesce in less than 1 Gyr because of the 
emission of gravitational waves (GWs): 

/1 1 \ 3 / 4 

acw "i)(^r) < V(e)1/4 ' (23) 

where Mi j6 = Mi/10 6 M© and, to 4th order in e, 

F(e) ee (1 - e 2)-V2 + g e 2 + | e 4^j (24 ) 

(Peters 1964). Assuming as before that three-body inter- 
actions begin operating at separation a , and combining 
equations (22) and (23), we get 

/ \ (l+7)/(12-4 7 ) 

^-10 3 ^(3-7)(-^) M-VV (e) -V 4 . 
«GW a \t + qj 

(25) 

liaf/acw < 1) then a MBHB would coalesce in less than 1 
Gyr after interaction with a bound stellar cusp. As shown 
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in the "coalescence diagram" of Figure 10, eccentric, mas- 
sive, very unequal-mass binaries embedded in steep stellar 
cusps are favored to reach coalescence. The steeper the 
cusp, the wider the portion of the Mi — q plane where 
coalescence can be reached within 1 Gyr. 
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Fig. 10. — Coalescence diagram. In each panel, the thick line 
delimits those MBHBs that can coalesce, because of GW emission, 
in less than 1 Gyr starting from separation aj (shaded area). As 
shown in the upper left panel, level curves to the left of such line 
are labelled according to the ratio a-f/acw reached by the binary 
at the end of the shrinking process (equation 25), while level curves 
to the right are labelled according to the the coalescence time at af . 



5. EVOLUTION OF THE STELLAR CUSP 

5.1. Cusp erosion 

The orbital decay of the pair occurs following the ejec- 
tion of a mass ~ 2 — 4 M2 (depending on cusp slope and 
binary eccentricity) and results in the progressive erosion 
of the stellar cusp. Eccentric binaries can shrink rapidly 
by scattering at pericenter a fewer number of deeply bound 
stars. More mass is expelled in the case of shallower cusps, 
because of the larger number of stars surrounding the bi- 
nary just outside do- Our hybrid approach allows us to 
compute the binary-driven evolution of the stellar density 
profile, once isotropy is assumed. At this stage, we do 
not perform a self-consistent treatment of the evolution of 
the anisotropy of the stellar cusp. Self-consistent integra- 
tion of the orbital decay yields dN e - } / da* , and the number 
of stars that remain bound to Mi as a function of a* is 
simply 

dN hd ^ d/V* __ dN^ 
da* da* da* 
We assume that the stellar cusp remains isotropic, i.e. 
(e„) = 0.667. The a* domain is then subdivided in inter- 
vals Aa*, and for each interval a number of bound stars 
oc (dNbd/da*)Aa* is generated with angular momentum 
distribution oc L*. We then compute the probability of 



finding a star at distance between r and r + dr from Mi 
(this is proportional to the time spent at such distance 
along its orbit), and reconstruct the stellar density pro- 
file p(r). Figure 11 shows the profile before and after bi- 
nary shrinking. The cusp is eroded between ~ O.Olao and 
<~ 2ao, depending on cq. For r < ao, a SIS is flattened to 
p oc r~ - 7 , while a r -15 cusp becomes p(r) oc r -0 ' 5 . Such 
results are independent of the mass ratio q. 

As three-body interactions and stellar ejections tend to 
circularize the orbits of ambient stars, one may wonder 
about the validity of the assumption that the cusp remains 
isotropic. We checked that even assuming that all stars 
still bound to the binary at the end of the process were set 
on circular orbits (i.e., maximum tangential anisotropy), 
then the cusp slope would be flattened at most by a fac- 
tor ~ 0.3 (for 7 = 2, and q = 1/729) with respect to 
the isotropic case. Differences in the cusp slope ^0.1 are 
common for initially shallower cusps and/or larger mass 
ratios. We conclude that, depending on 7 and on the final 
anisotropy of the bound stars, the slingshot mechanism 
creates a central core as flat as p(r) oc r -(°- 4 ^ - 8 ). 




Fig. 11. — The cusp stellar density profiles (in arbitrary units) 
before (thin lines) and after (thick lines) erosion by the shrinking 
binary. Solid line: q = 1/9; Long-dashed line: q = 1/27; Short- 
dashed line: q = 1/81; Dot-dashed line: q = 1/243; Dotted line: 
q = 1/729. The upper (lower) set of curves is for 7 = 2 (7 = 1.5). 
Fiducial r~ 0,7 and r~° power-laws are shown for reference. 



5.2. Distribution of kick velocities 

The distribution of stellar ejection velocities in units of 
the binary circular velocity at the initial separation a , 
V c fi, is shown in Figure 12. The distribution cannot be 
fit by a single or a broken power-law, as in the case of 
the kick velocities imparted at fixed binary given separa- 
tion, as the derivative of the distribution is a monotonic 
decreasing function of the ejection speed. While the peak 
of the distribution shifts toward smaller V/V c .o values as 
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q decreases (Vp Ca k oc i/q), the high velocity tail is inde- 
pendent of q. For initially eccentric binaries, the velocity 
of Mi at periastron is > V C: o and a significant number of 
HVSs can then be generated. Note that the speed distri- 
bution is only weakly dependent on 7 once V is expressed 
in units of V c .q. 

5.3. Ejected mass 

Integration of the kick velocity distribution gives the 
ejected mass M ej . This quantity is plotted in Figure 13 
as a function of q for different lower velocity thresholds. 
As discussed above, the total number of ejected stars is 
approximately oc q and is weakly dependent on eo and 7. 
Stellar ejection occurs in a burst lasting from few tens to 
several thousands binary orbital periods (see Fig. 14, note 
the different scale of the time axis in the four panels), i.e. 
from 10 5 to 10 7 yr. The highest velocity stars are delayed 
with respect to the bulk of ejections, as large kicks require 
close binary separations. The ejection rate peaks at earlier 
times in the case of more eccentric binaries, and is larger 
for steeper cusps. 
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Fig. 12. — Final velocity distribution of ejected stars as a func- 
tion of V/V c fi. Short-dashed lines: q = 1/9; Long-dashed lines: 
q = 1/81; Solid lines: q = 1/729. The assumed values of 7 and eo 
are listed in each panel. 

It must be pointed out that the quantity M e j plot- 
ted in Figure 13 accounts only for the mass ejected by 
energetic three-body interactions occurring in the final 
stage of the binary evolution. During the entire binary 
evolution, stars in galaxy nuclei are also displaced by 
the heating associated with dynamical friction - the cu- 
mulative effect of many weak encounters with distant 
stars. According to Merritt (2006), the total mass dis- 
placed by dynamical friction depends only weakly on 
binary mass ratio (for q in the range 0.025 < q < 
0.5), and has little or no effect on the innermost stel- 
lar cusp for q < 0.1 (see Fig. 5 in Merritt 2006). 



Fig. 13. — Mass ejected from the binary, scaled to the total bi- 
nary mass, as a function of q. Left panel: 7 = 2. Right panel: 
7 = 1.5. Thick lines: eo = 0.1. Thin lines: eo = 0.9. In each set 
of curves lines refers, from top to bottom, to all stars with kick 
velocities V > (the total mass expelled by the binary) and with 
log (V/ V c , ) > -1 + 0.5n for n = 0, 1, 2, 3. 
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Fig. 14. — Mass ejection rate for binaries with eo =0.1 and differ- 
ent mass ratios embedded in an SIS. Lines refers, from top to bot- 
tom, to all stars with kick velocities V > (the total mass expelled 
by the binary) and with log (V/V c ,o) > — 1 + 0.5n for n = 0, 1, 2, 3. 
Note that the time axis has a different scale in each panel. 
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While the total mass deficit accumulated during binary 
evolution appears to scale with the mass of the binary 
(i.e., to Mi for small q), the mass ejected after the en- 
ergetic encounters discussed in this paper, and therefore 
the number of expected HVSs, scale with the mass M 2 of 
the intruder. The number of stars ejected above a given 
velocity threshold (in units of V c .o) is also a weak function 
of eo and 7. 

5.4. Angular properties of HVSs 

The angular properties of HVSs show several peculiar 
features, qualitatively similar to those discussed in Pa- 
per I for the scattering of an unbound stellar population. 
The ejected stars are flattened in the binary orbital plane, 
and, in the case of eccentric binaries, are grouped into a 
"broad jet" aligned to the velocity of M 2 at periastron. 
Both anisotropies are more pronounced in stars undergo- 
ing a stronger interaction and receiving larger kicks. The 
anisotropy of the HVS population can be quantified by 
computing as a function of binary separation the angles 
(9 2 ) and (</>), where — tt/2 < 9 < ir/2 is the latitude of the 
star (the angle between the velocity vector V at infinity 
and the binary orbital plane) , and < <f> < 2ir is its longi- 
tude (the angle between the projection onto the binary or- 
bital plane of the velocity vector at infinity and the x-axis) . 
For a spherically-symmetric distribution, (9 2 ) ~ 0.47, and 
(0) = 7T. We find (9 2 ) ~ 0.35 - 0.40 for all HVSs. As 
a general trend, stars ejected above a given speed tend 
to become more isotropic as the binary shrinks, confirm- 
ing the analytical result of Levin (2006). This effect was 
already seen in three-body scattering experiments of un- 
bound stars (see Paper I for a detailed discussion), and 
is related to the fact that, as the pair decays, its circular 
velocity grows and even stars experiencing relatively weak 
encounters can attain large kick velocities. In terms of 
longitude, HVSs are ejected almost isotropically, though 
an azimuthal anisotropy become apparent with increasing 
kick velocities. The high- velocity tail of the distribution is 
formed by stars expelled after a close encounter with M 2 
near its periastron, in a broad jet with {(f>) ~ 3/2tt. As the 
binary potential is non-Kcplcrian, such broad jet will pro- 
cess during binary evolution on a timescale that depends 
upon q and 7 (Levin 2006). 

6. THE CASE OF SGR. A* 

Most of the results presented in the previous sections 
are scale invariant, i.e. they are independent on the ab- 
solute value of M and on the chosen normalization of the 
stellar cusp. The only underlying assumption is that the 
interaction between the MBHB and ambient stars begins 
when the stellar mass inside the binary orbit is equal to 
2M 2 . We are interested in scaling our results to the scat- 
tering of stars bound to Sgr A*, the massive black hole 
in the Galactic Center, by an inspiraling companion of in- 
termediate mass (Yu & Tremaine 2002; Sesana, Haardt & 
Madau 2007b). Let us first express equations (20), (21), 
and (22) in physical units: 



where uioo is measured in units of 100 kms . From ao, 
the time unit of our experiments is then 



r =0.43 pc (3-7)M li6 ar 2 , 



(27) 



Po = 1.96 x 10 6 M pc" 3 (3 - 7 r 2 Mja 6 W0 , (28) 

a = 0.43 pc (3-7)M 1 , 6( rr 2 g 1/(3 - 7) (l + 'z) (2 " 7)/(3 " 7) , 

(29) 



P Q = 1.3x10 s yr (3- 7 ) 3 / 2 M 1 1/ 6 4 q 



1 + q 



(2 7 -3)/(6-2 7 ) 



(30) 

The stellar density profile around the Galactic Center 
can be described as a double power-law, with outer slope 
~ —2 and inner slope ~ —1.5 (Genzel et al. 2003; Schodel 
et al. 2007). The massive black hole Sgr A* weighs 
~ 3.5 x 10 6 M Q (Schodel et al. 2002; Ghez et al. 2005). 
Using Mi, 6 = 3.5, ctioo = 1 and 7 = 1.5, from equations 
(27) and (28) we obtain r = 2.26 pc, and p = 7 x 10 4 M 
pc -3 , in good agreement with the most recent observations 
(Merritt 2006; Schodel et al. 2007). For the relevant values 
of q and 7, the typical timescale for orbital decay ranges 
between ~ 10 5 and <~ 10 7 yrs (see eq. 30). In the fol- 
lowing, we will consider two different mass ratios for the 
putative MBHB at the Galactic Center, q = 1/243 and 
q = 1/729, corresponding to an inspiraling IMBH of mass 
M 2 ~ 1.4 x 10 4 M Q and M 2 ~ 4.8 x 1O 3 M , respectively, 
as well as two different slopes for the initial inner stellar 
cusp, 7 = 1.5 and 7 = 1.75. We will also study the im- 
pact of binary initial eccentricity. The parameters of the 
different models considered are listed in Table 2. 
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Table 2 

Parameters of two different models for the stellar cusp at the 
Galactic Center. The quantities 7, ro, po, q, 0,0, Po, and V Cj o are, 
respectively, the cusp slope, the cusp characteristic radius, the 
density at ro, the binary mass ratio, the binary separation at 
which gravitational slingshots start, the binary orbital period at 
ap , and the binary circular velocity at ao . 

The dynamical evolution of a putative IMBH-Sgr A* bi- 
nary is displayed in Figure 15 for e = 0.1. The two upper 
panels show the time-changing semi-major axis a and ec- 
centricity e. The former shrinks to 10~ 2 — 10~ 3 pc in 2-15 
Myrs, depending on q and 7. Lighter IMBHs reach smaller 
separations on a longer timescale. Our results for the case 
7 = 1.75 can be directly compared to the numerical simu- 
lations of Baumgardt et al. (2006) and Matsubayashi et al. 
(2007), and are found to be in excellent agreement. The 
eccentricity increases rapidly in all cases to values ;> 0.8. 
We have checked that, when eo > 0.3, the binary eccentric- 
ity can reach values as large as e ^ 0.95. This is again in 
agreement with the results of Baumgardt et al. (2006) and 
Matsubayashi et al. (2007). In the figure we have marked 
with an horizontal dotted line the separation at which the 
binary can coalesce in 1 Gyr because of GW emission (eq. 
23). Only for a cusp as steep as 7 — 1.75 does the pair ac- 
tually reach such separation. One should note that, while 
in our hybrid model the eccentricity evolves smoothly, in 
a realistic situation it will undergo discontinuous "jumps" 
triggered by rare close encounters, which could induce ex- 
treme eccentricities and accelerate coalescence. The two 
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lower panels of Figure 15 show the resulting stellar den- 
sity profiles after binary erosion. Cusps arc flattened to 
p oc r~ ' 7 in the central fewxlCP 2 pc. The numerical sim- 
ulations of Baumgardt ct al. (2006) and Matsubayashi et 
al. (2007) produce somewhat shallower slopes, a discrep- 
ancy that may be associated with our assumption of an 
isotropic stellar cusp after the interaction. 
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Fig. 15. — The case of Sgr A*. Upper panels: time evolution of 
binary semi-major axis (solid lines, left axis scale) and eccentric- 
ity (dashed-line, right axis scale) for eo = 0.1. The dotted lines 
mark the separation at which binaries can coalesce because of GW 
emission in < 1 Gyr. Lower panels: evolution of the stellar density 
profile. In all panels, the thin lines are for q = 1/243, the thick lines 
for q = 1/729. 

The heating of the cusp results in the creation of a pop- 
ulation of HVSs. This is of particular interest since the 
discovery of the first HVS in the Milky Way (Brown ct 
al. 2005), and an IMBH inspiral onto Sgr A* is regarded 
as a possible source of hyper-velocity ejections (though 
the tidal break-up of close binaries by Sgr A* seems to 
be supported by observations and statistical studies, see 
e.g. Perets 2007, Sesana et al. 2007b). In Figure 16 we 
plot the stellar ejection rates as a function of time for dif- 
ferent models and for different velocities at the radius of 
influence of Sgr A*. Modelling the Milky Way potential 
as the sum of a luminous component (Miyamoto & Nagai 
1975) and a dark matter halo (Widrow & Dubinski 2005), 
we find an escape velocity from the Milky Way of ~ 840 
km/s at rj n f. This velocity threshold translates, in such 
a gravitational potential, into about 450 kms" 1 10 kpc 
away from Sgr A*. Stars with V > 300 kms' 1 at n n f do 
not leave the bulge, while stars with V > 600 kms -1 can 
reach 4 kpc away from the Galactic Center. Stars with 
V > 900 kins -1 and V > 1200 kins -1 are not bound to 
the Milky Way and, at a reference distance of 10 kpc, have 
still velocities of 600 kins -1 and 1000 kms -1 , respectively. 
Figure 16 shows that the ejection occurs in a relatively 
short burst lasting a few Myrs, with higher velocity stars 
being produced at increasingly later times. At peak, the 



ejection rate of HVSs with V > 900 kms^ 1 varies between 
5 x 10 -5 and 2 x 10~ 3 yr _1 , depending on q and 7. The 
rate is larger in steeper cusps, as stars are more centrally 
concentrated and are scattered when the binary separation 
is smaller and orbital velocity larger. 
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Fig. 16. — Ejection rate of HVSs for an IMBH falling onto Sgr 
A*. Solid lines: eo = 0.1; dashed lines: eo = 0.9. For each set of 
curves, lines from top to bottom correspond to a velocity threshold 
V > 300, 600, 900, and 1200 kms" 1 at the radius of influence of 
Sgr A*, respectively. 

The speed distribution of HVSs depends on the details 
of the model. In the range 300— 1000 kins" 1 (the velocity 
range of the HVSs observed by Brown ct al. 2006), and at a 
Galactocentric distance of 55 kpc (the average distance of 
the observed HVSs), the distribution can be approximated 
by a power-law f(V) oc V^ 1 - 5 , almost independent of q, eo, 
and 7. Assuming to* = 1 Mq, we find a predicted number 
of HVSs with V > 840 kms" 1 at r inf A^hvs - 525 (1290) 
for q = 1/729 and 7 = 1.5(1.75). This number roughly 
doubles for q = 1/243 and is fairly independent on e . Our 
results are consistent with Baumgardt et al. (2006), who 
estimate AWs ~ 1700 for M 2 = 10 4 M , and 7V H vs ~ 900 
for M 2 = 3 x 10 3 M Q . A peak of ejection occurs after 
^1 — 2 Myr. Levin (2006) finds comparable numbers of 
stars expelled, and a similarly peaked ejection rate. 

7. SUMMARY 

We have performed, for the first time, scattering exper- 
iments between a MBHB and stars drawn from a cusp 
bound to the primary hole. We have studied the dynam- 
ics of the pair and its orbital decay by three-body inter- 
actions, the impact of the gravitational slingshot on the 
stellar density profile, the properties of the ejected stellar 
population, and have scaled our results to the case of Sgr 
A*. Our results can be quickly summarized as follows: 

1 . The extraction of the cusp binding energy causes the 
binary to shrink by a larger factor compared to the 
scattering of unbound stars. The effect is more no- 
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ticeable in the case of small mass ratios q. 

2. The binary orbital eccentricity increases much more 
rapidly compared to the unbound case. The eccen- 
tricity growth is more pronounced in small mass- 
ratio binaries, and for shallower stellar cusps. 

3. The combined effects of enhanced orbital decay and 
eccentricity growth lead very unequal-mass binaries 
to the gravitational wave coalescence phase. The de- 
tailed fate of the pair depends on the absolute value 
of its mass. More massive binaries decay faster. 

4. The stellar cusp is eroded, and the total mass re- 
moved by strong three-body encounters is 2-to-4 
times the mass of the secondary hole. While the 
mass deficit caused by dynamical friction in a merger 
event scales with the binary mass (and involves 
mostly distant stars), the mass of stars ejected from 
the inner cusp by highly energetic interactions scales 
with M 2 . Ejection occurs in a "burst" lasting from 
few tenths to several thousands binary orbital peri- 
ods, depending upon q. 



5. Scaled to the scattering of stars bound to Sgr A* by 
an inspiralling IMBH, our results imply the forma- 
tion of a core of 0.1 pc in 1-10 Myrs, as well as the 
ejection of 500-2500 HVSs moving with speeds suf- 
ficient to escape the gravitational field of the Milky 
Way. In Sesana et al. (2007b) we have used the 
Brown et al. (2007) sample of unbound and bound 
HVSs together with numerical simulations of the 
propagation of HVSs in the Milky Way halo to con- 
strain this ejection mechanisms, and shown that it 
appears to produce a spectrum of ejection veloci- 
ties that is too flat compared to the observations. 
Future astrometric (as, e.g., GAIA) and deep wide- 
field (as, e.g., LSST) surveys should unambiguously 
identify the ejection mechanism of HVSs, and probe 
the Milky Way potential on scales as large as 200 
kpc (Gncdin ct al. 2005; Yu & Madau 2007). 
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